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We report the results of Molecular Dynamics simulations of electron transfer activation parameters 
of plastocyanin metalloprotein involved as electron carrier in natural photosynthesis. We have 
discovered that slow, non-ergodic conformational fluctuations of the protein, coupled to hydrating 
water, result in a very broad distribution of donor-acceptor energy gaps far exceeding that observed 
for commonly studied inorganic and organic donor-acceptor complexes. The Stokes shift is not 
affected by these fluctuations and can be calculated from solvation models in terms of the response 
of the solvent dipolar polarization. The non-ergodic character of large-amplitude protein/water 
mobility breaks the strong link between the Stokes shift and reorganization energy characteristic 
of equilibrium (ergodic) theories of electron transfer. This mechanism might be responsible for low 
activation barriers in natural electron transfer proteins characterized by low reaction free energy. 



I. INTRODUCTION 

Redox proteins play diverse roles as electron carriers 
in biological energy chains. 1 Enzymatic activity often in- 
volves transferring electrons to carry chemical reactions, 2 
while metalloproteins deposited in mitochondrial mem- 
branes and photosynthetic units serve as redox sites with 
tuned redox potential to allow one-directional electron 
flow in electron transfer chains; 3 Plastocyanin (PC) from 
spinach is a single polypeptide chain of 99 residues form- 
ing a /3-sandwich, with a single copper ion coordinated 
by 2 sulfurs from cysteine and methionine and 2 nitro- 
gens from histidine residues (Figure [1]). The presence 
of the copper ion, which can change redox state, allows 
PC to function as a mobile electron carrier in the pho- 
tosynthetic apparatus of plants and bacteria. It accepts 
an electron from ferrocytochrome / and diffusionally car- 
ries it to another docking location where the electron is 
donated to the oxidized form of Photosystem l4 

This functionality is achieved through fast electron 
transfer reactions at docking locations with low driving 
force ~ 20 meV and electron tunneling distance > 10 
Aj^ The efficient turnover of the photosynthetic appara- 
tus demands fast rates at redox sites, faster than typical 
biological catalytic rates of 10 2 — 10 4 s^ 1 (ref d). Given 
the small driving force, this constraint limits the reor- 
ganization energy A to about 1 eV»^ The reorganization 
energy here is a sum of the solvent, A s , and internal, Aj, 
components, where A s generally incorporates the com- 
bined electrostatic effect of the protein and water. In the 
rest of the paper, we will separate the atoms with partial 
charges varying with the redox state as the redox site 
(Figure [J) , considering the rest of the protein and water 
as the thermal bath. Our main focus will, however, be on 
the interaction of the redox site with water and that is 
how we define the solvent reorganization energy A s sep- 
arating the interactions with the protein atomic charges 
into the protein reorganization energy A pro t (see below 
for more precise definition). 

Both Xi and A s are large in synthetic redox systems 
with a copper ion serving as the redox site because 




FIG. 1: Structure of plastocyanin and the illustration of 
the protein large-scale conformational motions displacing hy- 
drating water. The active site includes copper ion (green), 
2 histidines (blue), methionine (red), and cysteine (orange) 
residues. The arrows and transparent parts of the protein il- 
lustrate motions of the main chain loops (not from actual MD 
simulations) displacing water molecules. 



of large structural changes upon electron transfer and 
a strong electrostatic interaction of copper with polar 
solvents. Experimental measurements^ and quantum 
calculations^ of the internal reorganization energy are 
still inconclusive placing it between 0.1 eV— and 0.6- 
0.7 eV (ref [6j and references therein). In addition, re- 
cent numerical simulations of heme and copper proteins, 
have uniformly placed their solvent reorganization ener- 
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gies in the range of 0.5-1.0 e Vi2iLSi2ii2iil These calcu- 
lations give results somewhat higher than what follows 
from the experimental work on Ru-modified aeruginoza 
azurin which has shown that the activation barrier dis- 
appears at A ~ 0.6 — 0.8 eVrH Even if the internal re- 
organization energy is as low as 0.1 eV, the available 
data suggest that electronic transitions involving cop- 
per proteins are significantly constrained thermodynami- 
cally requiring a tight docking configuration 13 and strong 
electronic overlap within the donor-acceptor pair which 
can therefore be modulated by protein's conformations. 14 
One therefore wonders if there are possibly some mecha- 
nisms at play, which are not included in standard models 
of electron transfer^ but which might allow a greater 
tolerance in varying the parameters affecting the activa- 
tion barrier. Our simulations reported here in fact show 
that the combination of charged surface residues with the 
coupled protein/water dynamic s 16 ' 17 ! 18 leads to a lower 
activation barrier without requiring either a larger driv- 
ing force or a higher electronic overlap. PC is used here 
as a prototype of what may apply to many other proteins 
involved in electron transfer chains given the wide spread 
of the type of protein/ water fluctuations considered here 
among other proteins not necessarily involved in redox 
activity^ 

II. ENERGETICS OF ELECTRON-TRANSFER 
ACTIVATION 

Electron transfer reactions are driven by thermal fluc- 
tuations of the nuclear degrees of freedom interact- 
ing with the electronic states of the donor and accep- 
tor. Electrostatic interactions between the atomic par- 
tial charges of the redox site with the partial charges or 
multipoles of the thermal bath usually follow the rules 
of the linear response approximatio n 10 ' 20 ' 21 embodied in 
Marcus theory of electron transfer.— The activation bar- 
rier is calculated in this picture from the crossing of two 
parabolic free energy surfaces Gi{X) depending on the 
donor- acceptor energy gap X. The use of equilibrium 
statistical mechanics to calculate Gi{X) results in several 
fundamental equations between the (spectroscopically 22 ) 
observable parameters of the model. The difference be- 
tween equilibrium vertical energy gaps (AX, Stokes shift) 
is equal to twice the reorganization energy X p and is also 
related to the variance of the energy gap a 2 — ((SX) 2 ) 
(spectral width): 

AX = X 01 - X 02 = 2X P = a 2 Jk B T (1) 

In eq [TJ X p refers to the solvent reorganization energy 
which, in traditional theories, is associated with the sol- 
vent polarization field.— This new notation is used here 
to distinguish the traditional definition of the solvent re- 
organization energy from our results for the solvent re- 
organization energy X s discussed below, which includes a 
new component not identified in the previous studies. In 
addition to eq[TJ energy conservation within Boltzmann 



statistics requires a linear relations^ 3 - 

G 2 (X) = G 1 (X)+X (2) 

Many attempts^ - including those for redox 
proteins^ ' 10 ' 21 to test eqs [1] and [5] have given pos- 
itive results validating the picture of two crossing 
parabolas. In contrast, our results here report a 
breakdown of both relations by coupled protein-water 
fluctuations affecting the statistics of the donor-acceptor 
energy gap. 

Protein electron transfer adds the protein matrix as 
an additional thermal bath characterized by a spectrum 
of vibrational modes including fast molecular vibrations 
incorporated into the internal reorganization energy Xi 
and slow conformational modes affecting the electrostatic 
potential at the redox sitej 24 ' 25 ' 26 ' 27 ' 28 This complication 
requires describing the reaction kinetics in terms of a 
multidimensional reaction coordinate space. The physics 
of the classical motions in the system is captured by a 
two-dimensional paraboloid energy surfac o 29 ' 30 ' 31 ' 32 ' 33 as 
a function of classical solvent, X, and effective vibra- 
tional, q, reaction coordinates (Figure H]). When both 
modes are fully equilibrated on the reaction time-scale, 
the reaction path Y = X + 75 is given as a linear com- 
bination of X and q with 7 representing the electron- 
phonon coupling. This reaction path dissects the two- 
dimensional space along the line connecting the minima 
of two paraboloids. The energetic separation between 
the minima defines the full Stokes shift AY related to 
the overall thermal dissipation of the energy of electronic 
excitations by the thermal bath. Extending eq[T]to equi- 
librium statistics in multidimensional coordinate space, 
one can obtain Ay in terms of the reorganization ener- 
gies: 

AY = 2(X P + A,) (3) 

However, when one of the modes is slow, the reaction 
path deflects from the line connecting the two minima 
and follows the fast reaction coordinate. The final state 
of the reaction then falls on the X-axis and is denoted 
by Xq2 in Figure [2] This picture, in which the solvent 
is a fast mode and the solute conformational mobility is 
a slow coordinate, was first considered by Agmon and 
Hopfieldi^ The problem of two-dimensional dynamics 
was later formalized by Sumi and Marcus who focused, 
in contrast, on the opposite case of fast intramolecular 
vibrations^ - 

The fully equilibrated path along the coordinate Y 
represents the lowest potential barrier between the two 
equilibrium points. If conformational equilibrium is not 
achieved on the reaction timescale tet = ^st (^et is 
the electron transfer rate), the reaction follows the path 
along X with the transition state marked by the cross 
on the X-axis (Figure [2]). However, if the reaction path 
deviates from the straight line due to stochastic confor- 
mational motions of the protein, it potentially can pass 
through a lower transition state marked by the cross on 



3 




FIG. 2: Electron transfer activation in two-coordinate space 
including solvent coordinate X and classical conformational 
coordinate q (multiplied by the factor a of electron-phonon 
interaction). The line connecting the two minima of the two- 
dimensional paraboloids corresponds to the reaction path Y 
for the fully thermalized fluctuations of both X and q co- 
ordinates. The total Stokes shift AY is then the energetic 
distance between the minima. Slow non-ergodic fluctuations 
of q shift the reaction path from the straight line connecting 
two equilibrium points Xoi to the wiggled line. The transition 
state then shifts from the cross point on the X-axis to a new 
point on the wiggled line. 



the wiggled line. The result of this is the breakdown of 
the link between the Stokes shift along the coordinate 
X, given as AX, and the effective curvature of the free 
energy surface determined by the variance of the energy 
gap fluctuations a 2 — ((5X) 2 ) (eq[TJ|. 

The modulation of the donor-acceptor energy gap by 
protein motions can be modeled by stochastic noise in 
contrast to equilibrium distribution resulting in eq [3J 
This effect is accounted for by adding average over con- 
formational fluctuations (subscript "q") to the Gaussian 
distribution along the solvent reaction coordinate 
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(4) 



Here, the dependence on q comes to the vertical en- 



ergy gap X 0i {q) 



i(q) * AEq (the asterisk refers 



to both the scalar product and space integration). This 
gap is formed by equilibrium solvent (nuclear) polariza- 
tion P eq i (g,r) in response to all partial charges of the 
protein (1439 atomic charges for PC) and the difference 
electric field AEo(r) created by the difference charges 
Azj of the redox site (Table SI, j runs over the atoms 
of the redox site). On the contrary, the polarization re- 
organization energy X p is calculated as the solvation free 
energy of Azj charges only, and is affected by positions 
of only a few atoms of the active site (4 in our simula- 
tions, see Sec. IIIip . Therefore, one can expect that it 
is the vertical gap that is predominantly modulated by 
protein motions while X p is mostly insensitive to such 
fluctuations (see below). 

Assuming Gaussian statistics of Sq and a linear expan- 
sion of X 0l (q) in 8q (X oi (q) ~ X oi + FSq), one gets: 



Gi(X) — Got + 



(X-X, 



0/ ) 



4(A p + A 9 (T obs )) 



(5) 



The new reorganization energy A g (T b s ) in principle car- 
ries the dependence on the redox state (i — 1,2), which 
requires non-parabolic energy surfaces^ and is not con- 
sidered here. 

The reorganization energy A g (r b s ) carries the depen- 
dence on the observation time r b s in order to stress on its 
non-ergodic characte r 32 ! 34 contrasting with equilibrium 
averages referring to r b s — > oo. The necessity to consider 
nonergodic activation parameters arises from the wide 
spectrum of relaxation times typical of proteins. A g (T bs) 
arises from the protein motions fast enough to produce 
energy gap fluctuations on the time frame r b s used to 
collect the averages. It can be obtained as the frequency 
integral of the autocorrelation function C q (uj) = (|<5(? w | 2 ) 
of 5q u with the low-frequency cutoff reflecting the final 
observation time 3 ^ 



A g ((T obs ) = {F 2 /k B T) / C,(w)dw 



(6) 



It turns into equilibrium reorganization energy X q — 
F 2 /(2k) (k is an effective force constant of harmonic con- 
formational motions) in the limit r b s — > oo considered 
by statistical mechanics. 

The conformationally-induced variance of the donor- 
acceptor energy gap 



CTq(7kin) = 2k B T kin X q (Ti 



kin) 



(7) 



is in principle accessible experimentally from hetero- 
geneous electron-transfer kinetics measured on pro- 
teins cryogenically quenched in their conformational 
substates. 2 - Here, the temperature of kinetic arrest Tk; n 
is estimated by requiring that the quenching rate Q — 
dT/dt and the temperature derivative of the conforma- 
tional relaxation time r q produce unity in their prod- 
uct: Q x (dr q /dT) = This approach, however, elimi- 
nates the hydration dynamics facilitating conformational 
changes (see below) . One can therefore expect that such 
experiments will inevitably underestimate a 2 observed at 
high temperatures. 

The arguments presented in this section are not meant 
to give an accurate theoretical description of the com- 
plex non-ergodic kinetics of electron transfer influenced 
by protein/water dynamics. They are more intended to 
set up a framework to understand the results of MD sim- 
ulations which provide a more detailed picture of the 
nuclear modes involved in the modulation of the donor- 
acceptor energy gap. 



III. COMPUTATIONAL METHODS 

A. MD Simulations 

Amber was used for all MD simulations. The 

initial configuration of PC was created using X-ray crys- 
tal structure at 1.7 A resolution (PDB: lag6 37 ). The 
system was heated in a NVT ensemble for 30 ps from 
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K to the desired temperature followed by volume ex- 
pansion in a 1 ns NPT run. After density equilibration, 
NVT production runs lasting from 15 ns (at 310 K) to 
18 ns (at 285 K) were made, of which 10 ns at the end 
of each trajectory were used to calculate the averages. 
The length of simulations was determined by monitoring 
the convergence of the solvent reorganization energy A s , 
which is the slowest-converging energetic parameter cal- 
culated here. The timestep for all MD simulations was 2 
fs, and SHAKE was employed to constrain bonds to hy- 
drogen atoms. Constant pressure and temperature sim- 
ulations employed Berendsen barostat and thermostat, 
respectively. 38 The long-range electrostatic interactions 
were handled using a smooth particle mesh Ewald sum- 
mation with a 9 A limit in the direct space sum. The 
total charge for the protein was —9.0 for the reduced 
(Red) state and —8.0 for the oxidized (Ox) state. Each 
state was neutralized with the corresponding number of 
sodium ions and TIP3P model was used for water ^ 

Three atomic charging schemes were utilized to 
parametrize PC's redox site (Table SI). For the first 
parameter set, a chemically fake charging scheme was 
employed that uses typical Amber force field (FF03 40 ) 
for all standard amino acid residues, but assigns an in- 
teger charge to the copper center in the reduced and 
oxidized states (Ql). Second, a more accurate charg- 
ing scheme was based upon experimental spin densi- 
ties from Solomon's group for the copper and copper 
ligands^i Finally, a third charge distribution is com- 
pletely parametrized at the DFT level for the charges 
and force constants of the copper and ligand atoms and 
consistent with the Amber force field (Q3).~ In addi- 
tion, Amber FF03 parametrization 4 ^ was applied to all 
non- ligand residues (Q2). There were various numbers of 
TIP3P water molecules for each of the charge distribu- 
tions: 5,874 (Ql), 5,886 (Q2), and 4,628 (Q3). 

We ran separate simulations (ca. 5 ns) for each charg- 
ing scheme to find that the results are not strongly af- 
fected by the choice of atomic charges. This was also 
noticed in some other recent simulations^^ We have 
therefore implemented charge scheme Q2 in all simula- 
tions reported here since it presents a good balance be- 
tween being simple and realistic. 

Amber force field 4 ^ was also used for the ground state 
tryptophan. Charges in the L a excited state were taken 
from the literature;* 3 - This charge set was chosen because 
it gives a good agreement with ab initio calculations of 
the indole ring4^ NVT simulations of tryptophan were 
carried out for a total of 3 ns in a simulation box contain- 
ing 420 water molecules after 1 ns density equilibration 
using NPT protocol with a Berendsen barostat^ 8 - The 
Stokes shift correlation function 4 ^ simulated for trypto- 
phan was in excellent agreement with both the experi- 
mental data 4 ^ and previous computer simulations. 47 This 
model simulation was used as a testing tool for our anal- 
ysis of the Stokes shift dynamics of PC. 



B. Calculations of the solvation thermodynamics 

Calculations of the solvent reorganization energy and 
the solvent part of the reaction free energy gap were car- 
ried out by two methods: (i) non-local response function 
theory (NRFT) i 49 and (ii) dielectric continuum approx- 
imation implemented in the DelPhi program suited Di- 
electric constant of TIP3P water (e s = 97. 5^ 9 -) was used 
for the solvent continuum and e s = 1 for the protein. 
This latter choice was driven by our desire to compare 
continuum and microscopic calculations of solvation ther- 
modynamics since the latter does not assume any polar- 
ization of the protein. A full account of the algorithm in 
application to the calculation of the redox entropy of PC 
will be published elsewhere^ 1 - In addition, since TIP3P 
water is non-polarizable eoo = 1.0 was used for the high- 
frequency dielectric constant in the reorganization energy 
calculations. 

In short, the NRFT calculation scheme employs the 
linear response approximation to replace the solvation 
chemical potential with the variance of the solute-solvent 
interaction potential Vo S ' 52 ' 53 

-fxos = (2k B T)- 1 {{SV 0s ) 2 } (8) 

The subscript "0" in the ensemble average (■ ■ • )o refers 
to the fact that, in the linear response approximation, 
the spectrum of electrostatic fluctuations of the solvent 
is not perturbed by the electrostatic solute-solvent inter- 
actions. Therefore, the variance in eq[5]is calculated for 
a fictitious system composed of water solvent and the re- 
pulsive core of the solute with all solute charges turned 
off. This approximation is known to work well for dense 
polar solventS ) 49 ' 53 i 54 and the main problem of the theory 
development is how to calculate the response function of 
the polar solvent in the presence of the solute which ex- 
pels the dipolar polarization field from its volumei 48 i 55 
This problem can be solved by applying the Gaussian 
solvation modelSS resulting in the linear response func- 
tion (2-rank tensor) x[Xs, ^o] functionally depending on 
the self-correlation function of the dipolar fluctuations of 
the solvent Xs(k) an d the shape of the solute occupying 
volume Qo- Once this problem is solved, the solvation 
chemical potential is calculated as a 3D, inverted-space 
integral convoluting the electric field of the solute Eo(k) 
with the response function: 

-Ha. = -E *x[x s ,^o] *E5 (9) 

Here, the asterisk refer to both the k-integration and 
tensor contraction and Eq is the complex conjugate of 
E . 

For polar liquids, the function Xs(k) splits into projec- 
tions longitudinal (parallel) and transverse (perpendicu- 
lar) to the wave-vector k. Each component is then rep- 
resented by the corresponding structure factor which is a 
function of the magnitude of k onl y 48 ' 49 These structure 
factors were obtained in this work from MD simulations 
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FIG. 3: Free energy landscape for conformational transitions 
(a) and a- relaxation (b) of the protein. G(Xo2,q) shows the 
cross-section of the two-dimensional free-energy surface taken 
at the equilibrium final value of the solvent reaction coordi- 
nate Xo2- Each conformational state along the reaction co- 
ordinate q contains a large number of substates separated by 
smaller free energy barriers (b). Transitions between these 
states are responsible for a-relaxation of the protein slaved to 
water dynamics. Each of a-substates can be separated into 
/3-substates (not shown) responsible for /3-relaxation of the 
protein and the hydration shell. 

of TIP3P water— at different temperatures (see refs l49l 
and [5TI for more details). With this input, the NRFT 
calculation is performed by grid summation (in k-space) 
of the solvent response function with the solute electric 
field. This latter was calculated numerically by using 
Fast Fourier Transform on the real-space lattice of 512 3 
points with a grid spacing of 0.42 A. In case of reorga- 
nization energy calculations, Eo(k) in eq[9]is obtained 
by taking only Azj charges of the redox site, thus pro- 
ducing the field AEo(k). In contrast, the solvent compo- 
nent of the free energy gap of electron transfer, AG S , is 
calculated with the complete charge distribution of the 
protein. 



IV. RESULTS 

The conformational dynamics of proteins are very 
disperse^ 7 including several time-scales that can poten- 
tially affect the energetics of electron transfer. Global 
conformational changes of the protein, which often occur 
on the time-scale of microseconds,— make the slowest 
time-scale. These transitions occur between minima of 
the free energy landscape separated by highest barriers. 
In terms of the two-dimensional coordinate space used 
in Figure [2j this motion sets up a transition from q = 
to q = Aq along the generalized protein coordinate q. 



In Figure [3] we show the cross-section of the free energy 
surface, F(Xo2,q), at the final state along the solvent 
polarization coordinate X = X a2 - The activation bar- 
rier separating the states q — and q — Aq is very high 
to be observed on the time-scale of our MD simulations. 
The topology of the free energy landscape^ is, however, 
more complex than that is sketched in Figure [3^- Each 
of the conformational states, q = and q = Aq, contains 
a large number of conformational substates separated by 
lower barriers^ (Figure [3Jd) . Transitions between these 
substates represent a-relaxation of the protein with many 
features analogous to a-relaxation of structural glasses. 60 
These protein dynamics are "slaved" to the solvent in 
a sense that the temperature dependence of the corre- 
sponding relaxation time follows that of water.— One 
of the consequences of this slaving is that the long- 
known dynamical transition of protein atomic displace- 
ments above the linear regime at T tr ~ 200-250 K±1£L£2. 
can be traced back to the fragile-to-strong dynamic tran- 
sition of hydrating waterj 18 ' 61 An alternative explana- 
tion suggests the merger of the fast f3- with slow non- 
observable a-relaxation at the transition temperature.'SS 
This fast /3-relaxation of the protein and hydrating water 
can be visualized as transitions between low-barrier sub- 
states within each landscape basin of the a-relaxation 
processes (not shown in Figure [3Jd) ■ Fast fluctuations 
between these substates involve amino-acid side chains 
and hydrogen-bond network at the protein surfac o 16 ! 63 
as well as protein vibrations which are not affected by 
the dynamic transition at T = T tr . /3— relaxation of the 
protein is strongly dominated by /3-relaxation of the hy- 
drating water— involving translational motions of water 
molecules in and out of the first hydration layer .A 6 - 

This scenario is consistent with the dynamics of the 
donor-acceptor energy gap observed along the MD simu- 
lation trajectory. A large-amplitude, redox-induced con- 
formational transition, if it exists^ is too slow to occur 
on the observation timescale r b s determined, in the com- 
puter experiment, by the length of the simulation trajec- 
tory. However, both a- and /3-relaxation of the protein 
and water (about 40% of the overall protein relaxation on 
the 10 ns time-scale^) are clearly seen in the X(t) tra- 
jectory (Figure|4]). The slower a-relaxation component is 
represented by large-amplitude oscillations superimposed 
onto fast /3-fluctuations of the solvent dipolcs. The slow 
a-relaxation is not typically seen in small rigid solutes ex- 
emplified by the trajectory of tryptophan (inset in Figure 
H]) where only fast /3-fluctuations are present. The time- 
scale of a- fluctuations (ca. 1 ns) suggests their origin in 
the motion of polar side groups^ which show jumps in 
their dihedral angles on the same time-scaled 



A. Interactions with water and protein 

The overall donor-acceptor energy gap AE is a sum of 
a gas-phase component, mainly affected by the ligands 
coordinating copper, and electrostatic interactions with 
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Simulation Time/ps 

FIG. 4: Trajectory of the solvent component of the donor- 
acceptor energy gap of PC in TIP3P water at 310 K. The 
upper trajectory shows unrestricted protein/water dynamics 
and the lower curve refers to protein atomic displacements 
frozen by applying positional harmonic restrains (8.0 kcal 
moP 1 A~ 2 ). Inset shows the same property for photoexci- 
tation of tryptophan in TIP3P water. The gray regions in 
the upper trajectory indicate segments of the trajectory 100 
ps long. The lower trajectory was shifted down by 1.5 eV for 
better visibility. 

the protein matrix, AE pTO t, and with the solvent (water), 
X. Their sum makes the total reaction coordinate 

Y = X + AE pmt (10) 

discussed in Sec. Mil In eq [TUl A_E prot is connected to 
characteristic conformational coordinate of the protein q 
by electron phonon coupling 7: .Eprot = 79- Table [J lists 
the statistics of Y fluctuations obtained from 10 ns of 
simulation data. We report the average donor-acceptor 
energy gap from the interaction of Azj charges (Table 
SI) of the redox site (Red and Ox states) with water, 
(X), and protein, (A_E prot ), as well as their sum, (Y). 
In addition, Table U gives the variance of Y which can 
be split into two self-correlation functions and one mixed 
protein- water component: 

a 2 = ((SY) 2 ) = a 2 + <7 2 pmt + a s , pmt (11) 

The splitting of the total variance into the reorganiza- 
tion energies A s = a 2 /(2k B T) and A prot = cr^ rot /(2fc B T) 
introduced in Sec.|T]neglects the cross-correlation compo- 
nent c^prot which, in our simulations, amounts to 5-18% 
of the total variance. This is the error bar for the separa- 
tion of protein and solvent nuclear fluctuations into two 
separate stochastic processes. 

B. Solvent reorganization energy 

Reorganization energy A g (r b s ) in eq [5] originates from 
fluctuations of the solvent dipolar polarization induced 
by coupled protein-solvent dynamics. This reorganiza- 
tion energy thus represents a new solvent mode affecting 
the donor-acceptor energy gap absent in the traditional 
theories of electron transfer operating in terms of sepa- 
rate vibrational (Aj) and polarization (A p ) nuclear modes. 
The enhancement of the energy gap variance by this new 
mode is very significant: the variance of X changes from 



a 2 = 2k B TX p (\ p = AX/2 ~ 0.45 - 0.65 eV), com- 
parable to other simulations, 7,21 to a much higher value 
a 2 = 2k B TX s characterized by the solvent reorganization 
energy A s = A p + A g ~ 5 eV (Table HI]) . This gigan- 
tic value of the reorganization energy far exceeds what is 
typically observed for electron transfer reactions between 
small hydrated ions.— 

The variance of the donor-acceptor energy gap a q ~ 
0.5 eV produced by conformational flexibility in our sim- 
ulations is significantly higher than experimentally re- 
ported o q ~ 0.05 eV from charge recombination in bac- 
terial reaction centers trapped by cooling (eq[7])in their 
conformational substates3i As noted above, this lower 
variance of energy gaps in reaction centers is expected 
since water dynamics significantly contributing to fluctu- 
ations of the donor-acceptor energy gap is also quenched 
by cooling. In addition, the hydrophobic environment of 
cofactors located in the membrane protein complex and 
the low temperature of the kinetic arrest 7ki n ~ 175 K 
(eq[7]) both contribute to the lower <j q . 

We need to emphasize that the reorganization energies 
considered here refer to the change in the charge distribu- 
tion of PC only. A donor-acceptor complex composed of 
two proteins (photosynthetic electron transfer) will also 
include a change of charges on the partner (heme) pro- 
tein. The interaction of this other set of charges located 
within a membrane protein with water is expected to be 
weaker than for hydrated PC. Therefore, there should be 
only minor Coulomb correction to the reorganization en- 
ergy. In addition, some reduction of the reorganization 
energy will arise from electronic polarizability of water 
not included in TIP3P parametrization (our calculations 
using non-local solvent respons o 48 ' 49 show a reduction 
from 0.74 eV in TIP3P water to 0.40 eV in ambient wa- 
ter). Nevertheless, the gigantic magnitude of A p + X q 
compared to the commonly considered X p calls for atten- 
tion to the effects of coupled protein/ water dynamics on 
electron transfer. 

The solvent effect on the electron transfer thermody- 
namics is dominated by water molecules closest to the ac- 
tive site. Protein flexibility significantly modulates this 
first solvation shell producing fluctuations of the closest 
Cu-0 distance around the average of 6.64 A, the largest 
fluctuation amplitude of ~ 2 A, and the standard devi- 
ation of 0.6 A (Figure El). We also note that the Cu-0 
pair distribution function (Figure S2) does not change 
with changing the redox state of plastocyanin, in contrast 
to observations reported for heme proteins i 7 ' 10 ! 11 ' 21 With 
such large-amplitude fluctuations, water is effectively fur- 
ther apart from the protein surface than the distance of 
the closest approach. The solvent part of the reaction 
free energy is then nearly 1.5 times smaller in magnitude 
than the value calculated from our solvation mode l 48 i 49 
assuming the closest approach (Table HT|) . In contrast, 
the calculated reorganization energy A p is in good agree- 
ment with simulations, which supports our assumption, 
used to derive eq\5\ that conformational fluctuations do 
not affect this parameter. 
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TABLE I: Averages (eV) and variances (eV 2 ) of the interaction energy between the Azj charges of the active site and the 
solvent (s) and protein (prot) for PC in Red and Ox states feci Hip . (Y) and a 2 refer to the total average energy gap and 
corresponding variance, respectively. The data are collected from a 10 ns trajectory. 

T/K Redox state (X) (AE plot ) (Y) a 2 Iot o'j <j' a q proM 



310 Ox 2.484 -8.587 -6.103 0.085 0.290 0.275 -0.051 

Red 3.394 -8.690 -5.296 0.078 0.298 0.416 0.014 

285 Ox 2.460 -8.465 -6.006 0.082 0.268 0.287 -0.053 

Red 3.773 -8.964 -5.190 0.071 0.233 0.249 -0.028 



TABLE II: Reorganization parameters of PC (all energies are 
in eV). 

T(K) Ap ~\q AG S a A pro t 6 

"285 08P 4.8(Red) -3.2 O 

0.77(0.69,3.6) d 4.1(Ox) -4.7(-7.1, -9.6) d 

310 0.54 c 5.0(Red) -2.9 1.5 

0.74(0.54,3.6) d 4.5(Ox) -4.6(-7.1, -9.6) d 

"Solvent component of the redox free energy obtained from the 
simulation data as AG S = Gf cd - G° x = -((X)i + (X) 2 )/2. 

'Calculated from the variance <r^ rot of the electrostatic interaction 
of Azj charges of the redox site with the charges of the protein as 

Aprot = CT 2 rot /(2fc B T). 

c Calculated from the simulation data as ((X)± — (X)2)/2. 

^Theoretical calculations using atomic charges. vdW radii, and 
coordinates of the protein combined with microscopic, non-local 
response functions of water ^3/19 Dielectric continuum calculations 
(brackets) have been done with the solvent-accessible cavity (first 
number) and standard vdW cavity (second number). 



We note in passing that X p from continuum calcula- 
tions is very sensitive to the definition of the dielectric 
cavity. When van der Waals (vdW) radii of protein atoms 
arc used to determine the cavity, high-polarity dielectric 
is allowed in a narrow pocket near copper thus signifi- 
cantly increasing the free energy of solvation. When, in 
contrast, a water molecule is rolled on the vdW surface 
to determine the solvent-accessible cavity, the result of 
continuum calculations is comparable to both the NRFT 
and MD numbers. We will provide a more detailed dis- 
cussion of these results in a separate publication^ 

The large value of \ q raises the question of whether 
the new nuclear mode responsible for the energy gap 
variation should be attributed solely to conformational 
motions of the protein or to more complex collective dy- 
namics coupling solvent to protein fluctuations. The evi- 
dence existing in the literature advocates the latter view 
suggesting that both the a- and /3-relaxation of the pro- 
tein are strongly coupled to hydrating water. Our at- 
tempts to connect the slow modulations of the X{t) tra- 
jectory (Figure^]) to the vibrational density of states of 
the protein^ have not given positive results since the 
low-frequency vibrations seen in Figure [4] could not be 
resolved from the quasi-harmonic analysis^ (Figure SI) 
or from the intermediate scattering function^ We have 
also tried to freeze the protein motions through harmonic 
positional restraints on atomic translations, with the re- 
straint weight equal to 8.0 kcal mol -1 A~ 2 . These simu- 
lations (ca. 5 ns started at the end of the unrestrained tra- 
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FIG. 5: Trajectory of the closest distance between Cu of 
the plastocyanin active site (Figure [TJ and oxygen of wa- 
ter in the Red (gray) and Ox (black) states. The inset 
shows time autocorrelation functions of the Cu-O distance 
in Red and Ox states. The autocorrelation functions fit well 
to eq [13] with the set of fitting parameters {Ag,tg,te, 0}'- 
{0.31, 0.2, 12.3, 0.53} for Ox and {0.39, 0.2, 1.1, 0.46} for Red. 
Here, tg and te are in picoseconds. 

jectory) have resulted in the energy-gap variance a 2 di- 
minished by a factor of ~ 3 (lower trajectory in Figured]), 
but still not reaching the value a 2 from the Stokes shift. 
This observation supports the view that ct-fluctuations 
of the donor-acceptor gap are coupled to translational 
motions within the hydration layer at the protein surface 
and that these fluctuations cannot be separated from pro- 
tein's conformational dynamics. Nevertheless, the strong 
reduction of a 2 upon freezing of the protein still suggests 
that protein motions produce the largest energetic con- 
tribution to the reorganization energy X q . 



C. Protein dynamics 

If the large- amplitude protein/ water motions affecting 
the solvent polarization are overdamped^ they can be 
described by Debye relaxation with an effective relax- 
ation time T q . The use of the Debye relaxation function 
in eq [S] gives the following simple equation for the non- 
ergodic reorganization energy^ 

A,j(T b s ) = (2A g /7r)cot- 1 (r 9 /T ob s) (12) 

Protein dynamics coupled with dipolar solvent polariza- 
tion are dominated by very slow motions with the char- 
acteristic time r q of about 0.5-1 ns, as follows from the 
fit of A p + A 9 (r b s ) (eq [T2"]) to the simulation data (Fig- 
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FIG. 6: Solvent (water) reorganization energy of Ox (open 
circles) and Red (filled circles) states of PC at 285 K vs the 
observation time r b s denned as the length of trajectory over 
which the averages were calculated. The solid lines are fits of 
simulations to eq (|12[) with r 9 = 1 ns (Ox) and r q = 0.5 ns 
(Red). The hatched diamonds indicate the results from ref 
[ly . The inset shows the initial portion of the plot. 

ure[6]). The non-ergodic component A q (T b s ) was calcu- 
lated from the energy gap variance taken on observation 
windows r b s along the simulation trajectory (gray seg- 
ments of length 100 ps in Figure[4]). On short observation 
times, r Q bs < 100 ps, i.e. fast electron transfer reactions, 
the slow conformational modulation does not show up, 
and the reorganization energy from the width approaches 
that from the Stokes shift, thus restoring eq[TJ A similar 
behavior, including the magnitude of the corresponding 
reorganization energy, was observed in ref (hatched 
diamonds in Figure [6]). 

With such slow conformational modulation of the wa- 
ter polarization, each short segment of the long tra- 
jectory finds itself in a different configuration, a situ- 
ation akin to dynamical heterogeneity responsible for 
stretched-exponential relaxation of structural glasses^ 
This picture is indeed confirmed by the Stokes shift 
correlation function calculated on segments of trajec- 
tory of different length. Stokes shift correlation function 
C(t) = (X(t)X(0)) from a short segment has a typical 
biphasic form composed of a fast Gaussian decay followed 
by exponential relaxation 4 ^ for which the stretch expo- 
nent (3 in eq[l3]is equal to unity: 

C(t) = A G e- (t / TG ^ 2 + (1 - A G )e~^/ TE ^ (13) 

On the contrary, the Stokes shift correlation function 
calculated on a longer segment (1-2 ns) develops a 
stretched-exponential relaxation with the stretching ex- 
ponent P = 0.69 and relaxation time of about 150 ps 
(Figure S3). This long tail, which may require longer 
simulations to be fully resolved^ is caused by collective 
water displacements (Figure [1]) by slowly moving parts 
of a biopolymer i 47 i 69 

D. Free energy surfaces 

The picture of non-ergodic, glassy dynamics emerging 
from the static and time-resolved energy-gap statistics 



is consistent with the free energy surfaces Gi(X). They 
are very shallow on the long 10 ns trajectory becom- 
ing increasingly curved on a shorter observation window 
(Figure \T&). The full free energy surfaces, obtained by 
sampling the total interaction energy of the active site 
with both the protein and the solvent (reaction coordi- 
nate Y in Figured]), ar e uniformly shifted to the negative 
values of Y without significant change in the Stokes shift 
(Tableland Figure [TJa). The slow non-ergodic dynamics 
of protein conformations thus decouples the Stokes shift 
from the reorganization energy breaking eq[T]down. Ac- 
cordingly, the difference G^AT) — Gi(A) in the region of 
overlap of Ox and Red surfaces is still a linear function of 
X, but with the slope of 0.10, instead of the unitary slope 
predicted by Marcus theory and usually observed in fully 
equilibrated systems J£ This number is consistent with eq 
[5]which yields the slope of AX/ (AX + 2\ q ) ~ 0.13. 

The statistics of the donor-acceptor energy gap in- 
duced by protein fluctuations are also approximately 
Gaussian. The width of the distribution is given by 
the variance of the electrostatic interaction energy of 
Azj charges of the redox site with the protein atomic 
charges. The corresponding reorganization energy is ap- 
proximately 1.5 eV (Table HI)) . However, since the reac- 
tion path is expected to follow the fast solvent coordinate, 
our main focus here is on the free energies Gi(X). 

These results do not contradict experimental estimates 
of A p + A, ~ 0.6 — 0.8 eV for copper proteins obtained 
from the top of the energy gap law when the reaction 
barrier disappears^ 1 ^ From eq[51 the activationless tran- 
sition is achieved at Aoi = when the reaction free en- 
ergy AGo obeys the equation — AGo — A; = AX/2 = \ p . 
With the inner reorganization energy currently estimated 
as low as 0.1 eV, 7 X p — 0.54 — 0.81 eV from present 
simulations is consistent with experiment. On the other 
hand, reorganization energy A s entering the distribution 
width a s is about an order of magnitude higher. The 
breakdown of the link between the Stokes shift and the 
distribution width (eq [T|) must have significant impli- 
cations for the biological function of PC and probably 
of other electron carrier proteins. The large distribu- 
tion width leads to an extremely small activation bar- 
rier, AG act = (X p + AG )7(4A P + 4A 9 ) ~ 0.08 eV, and 
thus fast electron transfer, for reactions at the docking 
locations with the typically small reaction free energy 
AGo — ~20 meV4 The standard rate estimate 3 then 
gives 13 A for the donor-acceptor distance at which the 
threshold catalytic rate of 10 4 s _1 is achieved. 



V. CONCLUDING REMARKS 

The application of the ideas presented here to biologi- 
cal electron transfer requires the transition from the ob- 
servation time determined by the length of the simulation 
trajectory to the reaction kinetics. This is achieved by 
setting T bs = tet which, given the activation barrier is 
a function of T Q b s , leads to a self-consistent equation for 
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FIG. 7: Electron transfer free energy surfaces of PC at 310 K 
plotted against the solvent reaction coordinate (a) and against 
the total (water+protein) reaction coordinate Y (b). The 
dashed lines in A are fits of Gi(X) from the 10 ns trajectory 
to eq ([5]). The narrow curve in (a) has been obtained by cal- 
culating the distribution functions on 100 ps segments of the 
trajectory (gray segments in Figure[4]| and averaging them af- 
ter sliding to a common maximum. All curves are logarithms 
of normalized distribution functions along the corresponding 
reaction coordinate. 

the rate^i 

fc ET oc exp [-AG act (£: E T)/(fcBT)] (14) 

Equation [TJJ incorporating the notion of reaction non- 
ergodicity, offers a compelling picture of the hierarchy 
of electron transfer reactions in photosynthetic systems. 
Faster reactions with /get ^ 10 9 s effectively cut off slow 
conformational motions of the protein from their energet- 
ics (Figure O, and the standard picture^ of the Stokes 
shift and parabolas' curvature related by eq [1] applies. 
Realizing fast electron transfer then requires activation- 
less transition as observed in primary charge separation 
in photosynthetic reaction centers (Figure [5Ji) On the 
contrary, reactions in the sub-nanosecond range start to 
experience the effect of conformational modulation of the 
activation barrier and can in fact proceed efficiently even 
with a small driving force since the activation barrier is 
lowered by the growing width of the energy gap distribu- 
tion (Figure [5b). Therefore, when speed is at stake, nat- 
ural systems have to lose in redox potential in exchange 
for fast activationless transitions. When slower reactions, 
still faster than catalytic rates, can be afforded, losing re- 
duction potential is not the necessity, and reactions with 
low driving force can still be efficient. 

Strong coupling between dipolar polarization and pro- 
tein mobility advocated here is consistent with the 
long suggested connection between protein dynamics 




X 



FIG. 8: Energetics of fast electron transfer reactions losing in 
redox potential to achieve activationless transitions (a) and of 
slower (ns range) reactions (b) which are allowed to proceed 
with a small driving force. The reaction barrier is lowered 
in the latter case by non-ergodic conformational/water dy- 
namics transforming the dashed-line parabolas into solid-line 
parabolas (b). The vertical arrow in (b) shows the suppres- 
sion of the activation barrier by reorganization energy X q . 

and hydration ! 16 i 18 i 19 i 61 i 62 The new reorganization en- 
ergy discovered here is related to the rms displacement 
of the slow mode as X q cx ((5q) 2 )/T and is therefore ex- 
pected to show a sharp increase at T > T tr when ((Sq) 2 ) 
starts its nonlinear rise. This observation might provide 
a resolution of the long-standing puzzle of electron trans- 
fer kinetics in many plants and bacteria: Arrhcnius plots 
of electron transfer rates often show breaks in their slopes 
at temperatures consistent with T tr i 71 ' 72 ' 73 This feature 
might be linked to the rise of \ at the onset of confor- 
mational activity in proteinsi^i 

The existence of the solvent-slaved a-relaxation is 
presently traced back to strong coupling between po- 
lar/ionized surface residues and water. This type of 
dynamics is usually not observed in inorganic and or- 
ganic donor acceptor complexes used for ET reactions 
and is presently believed to be unique to natural poly- 
mers. Properties of synthetic polymers, in particular in 
respect to glassy dynamics^ have many common fea- 
tures with biopolymers and one hopes that phenomena 
analogous to ones observed here might be realized in 
flexible donor-acceptor architectures including branched 
polymers and dendrimeric structures. 
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